import numpy as np
# pip install matplotlib
import matplotlib.pyplot as plt

Y = np.zeros((6, 30000,), dtype=np.float32)  # 30000代表30s
R = np.zeros((6, 30000,), dtype=np.float32)  # 30000代表30s
Δt = 0.001  # 时间间隔为1ms
ωnArray = [0.2,2,8]
ξ = 0.2
for j in range(0, 3, 1):
    ωn = ωnArray[j]
    yk_2 = 0
    yk_1 = 0  # 历史数据
    rk = 0
    for i in range(0, 30000):
        if (i > 999):
            rk = 1  # 阶跃信号

        a = 2*yk_1 + 2*ξ*ωn*yk_1*Δt - yk_2 + rk*ωn*ωn*Δt*Δt  # 分子
        b = 1 + 2*ξ*ωn*Δt+ωn*ωn*Δt*Δt  # 分母
        yk = a/b
        Y[j][i] = yk
        R[j][i] = rk
        yk_2 = yk_1  # 历史数据
        yk_1 = yk  # 历史数据

t = [i*0.001 for i in range(0, 30000, 1)]  # 5000个数
plt.title("Simulation ξ=0.2")  # 括号当中输入标题的名称
for k in range(0, 3, 1):
    plt.plot(t, Y[k], label="ωn=%.1f" % ωnArray[k])
# plt.plot(t, R[0], label="R")
plt.legend(loc='best')  # 图列位置，可选best，center等
plt.show()
